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,— i 1 Abstract. A thin gaseous disk has often been investigated in the context of various 

O j phenomena in galaxies, which point to the existence of starburst rings and dense cir- 

cumnuclear molecular disks. The effect of self-gravity of the gas in the 2D disk can be 
important in confronting observations and numerical simulations in detail. For use in 
such applications, a new method for the calculation of the gravitational force of a 2D 
disk is presented. Instead of solving the complete potential function problem, we cal- 
culate the force in infinite planes in Cartesian and polar coordinates by a reproducing 
kernel method. Under the limitation of a 2D disk, we specifically represent the force as 
a double summation of a convolution of the surface density and a fundamental kernel 
. and employ a fast Fourier transform technique. In this method, the entire computa- 

tional complexity can be reduced from 0(N 2 x N 2 ) to 0(N 2 (log 2 N) 2 ), where N is the 
number of zones in one dimension. This approach does not require softening. The pro- 
posed method is similar to a spectral method, but without the necessity of imposing a 
periodic boundary condition. We further show this approach is of near second order 
accuracy for a smooth surface density in a Cartesian coordinate system. 
AMS subject classifications: 52B10, 65D18, 68U05, 68U07 
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1 Introduction 

The potential O of a given distribution of density p in IR 3 satisfies the Poisson equation, 

AO(x)=47rGp(x)=/(x), xGR 3 , (1.1) 



* Corresponding author. Email addresses: yenOmath . f ju.edu. tw (C.C. Yen) 

http: / / www.global-sci.com/ Global Science Preprint 



2 



where G is the gravitational constant and x = (x,y,z) is the position. Without loss of 
generality, we may assume that the gravitational constant G=l. Provided that the density 
profile has a continuous second derivative with respect to the spatial coordinates, the 
potential is smooth. In this situation, the numerical approach for solving the potential 
via (|1.1|) by the finite difference method is adopted. Artificial boundary conditions are 
imposed in the numerical approach for solving dl.ll) because the boundary condition is 

lim <&(*) = 0. (1.2) 

The Poisson equation is intrinsically 3-dimensional, and the calculation of the potential 
can be computationally prohibitive. A possible solution to reduce the computation time 
is to apply the multigrid method 03 El/ but the computational complexity is 0(N 3 ), 
where N is the number of zones in one dimension. 

1 

The solution of (|1.1|) can be represented in terms of the fundamental solution, — fC(x), 
where 

1 

KLix 



^/x 2 +y 2 +z 2 ' 



as 



<$>(x,y,z) = — J J J1C(x — x,y — y,z — z)f(x,y,z)dxdydz. (1.3) 

The above formula is preferable to (|1.1)> when the density is not smooth. The potential 
can be solved via the integral equation in dl.31) . Spectral methods are a common method 
of choice and a review article has recently been written by Shen and Wang [8], describing 
work on the analysis and application of these methods in unbounded domains. The 
difficulties encountered in the numerical approach for solving dl.ll) or dl.3l) are related to 
the extent of the domain R 3 and the density which can be singular. 
In this paper, we consider the density represented by 

p{x)=a(x,y)5{z), (1.4) 

where o~(x,y) is so-called surface density equal to 

cr(x,y) = J p(x)dz. (1.5) 

We restrict our attention to calculating the forces directly for the surface density of com- 
pact supports. 

For an infinitesimally thin gaseous disk, the multigrid method, which is intrinsically 
suited for 3D problems, cannot be reduced for the two dimensional problem we consider 
in this paper. The spectral method using Fourier basis functions on a two dimensional 
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space artificially imposes the assumption of periodic boundary conditions. This is not re- 
alistic for the long range gravitational force calculations. A direct method without the pe- 
riodic assumption requires a softening parameter technqiue, but the accuracy is reduced 
simultaneously. A method is proposed which is of linear complexity without artificial 
boundary conditions, and near second order accuracy. 

This paper is organized as follows. The framework and assumption are presented 
in Section [2] Sections |3] and H] describe the numerical methods for Cartesian and polar 
coordinates, respectively. Section [5] demonstrates the order of accuracy of the proposed 
methods as verified by a family of finite disks (e.g., D2 disk; [71]) and a disk of a pair of 
spirals. A comparison with several existing methods is also presented in that section. 
Finally, the discussion and conclusion are given in section [6] 

2 Framework and assumption 

The evolution of a thin disk is of fundamental interest in astrophysics and the effect of 
the self-gravity of gas therein may be important in modeling observed phenomena in 
detail. This paper presents a numerical method for solving the self-gravitating forces 
in Cartesian and polar coordinates, which can be used in modeling infinitesimally thin 
disks in galaxies and protostellar systems ITOI . 

The self-gravitating force can be determined by taking derivatives of the potential 
function which satisfies the Poisson equation in dl.lfr . However, the calculation of the 
potential dl.lt is on an unbounded domain and the solution in a finite region requires 
the imposition of artificial boundary conditions. The solution of Poisson's equation with 
variable coefficients and Dirichlet boundary conditions on a two dimensional irregular 
domain is one of second order |2|. 

Let us confine our attention to the density in an infinitesimally thin disk as defined in 
dl.4|) and dl.5|) . Here, we focus on the self -gravitating force computation. The approach 
presented in this paper is to directly calculate the self -gravitating force by expressing the 
potential function as a type of a convolution of the surface density and the fundamental 
kernel and taking the derivative of the potential function. This approach is similar to 
the spectral method, but less restrictive. Trigonometric bases functions and the artificial 
periodic boundary conditions are used for the spectral method, but are not required in 
the proposed approach here. 

A uniform grid discretization in Cartesian coordinates and a linear approximation of 
the surface density on each cell are used to reduce the computational time and increase 
the accuracy of the numerical solution, respectively. Similarly, for polar coordinates, a 
logarithmic grid discretization is used instead of a uniform grid discretization. Based on 
the discretization and approximation, the self-gravitating force is written as a convolu- 
tion form of double summations. It is known that the calculation of convolution form can 
be accelerated by the use of a fast Fourier transform (FFT), see Appendix B. Employing 
the FFT, the computational complexity is reduced from 0(N 4 ) to 0((Nlog 2 N) 2 ), where 
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N is the number of zones in one direction. The linear approximation also leads to an order 
of convergence that is near second order 0(h 2 ), where the size of a zone h = 0(l/N). 

3 Self -gravitating force calculation in Cartesian coordinates 

In this section, we describe the method in detail. The potential function O of dl.lt can be 
expressed as 

<5(x,y,z) = — J J J K,(x — x,y — y r z — z)p(x,y,z)dxdydz, 
1 

where /C(x,y,z) = — =. By (11.4ft . the forces on the disk in the x-direction and the 

y / X 2 + t/ 2 +Z 2 
y-direction become 

— O(x,y,0) = / / —K(x-x,y-y,Q)a(x,y)dxdy (3.1) 

and 

3 f f d 

—<P(x,y,0) = J J — K{x-x,y-y,0)a{x,y)dxdy. (3.2) 

We calculate d3.lt and (|3.2D by a numerical approach. Here, we focus on the derivation of 
the force calculation in the x-direction. The force in the y-direction is obtained in a similar 
manner (see Appendix A). 

Since the support of the surface density is compact, contained in a domain D = 
[— M,M] x [— M,M] for some number M>0, we discretize the region uniformly as follows. 
Given a positive integer N, we define Ax = 2M/N, Ay = Ax, x !+ i/2 = —M+iAx, y/+i/2 = 
—M+jAy, where i,j=0,...,N. We further define the center of the cell D,- / = [x,_ 1 / 2/ x i+1 / 2 ] x 
[1//-1/2/I//+1/2] as Xi = (x,_i /2 + x i+ i/ 2 )/2 and y } = (y ; -_ a /2+yy + i/ 2 )/2 / where i,j = l,...,N. 
Hence, the domain D is discretized into the N 2 cells. 

The forces in the x-direction and the y-direction at the center of cells are 

F^. = Ao(x ! ,y ; ,0), and Ff. = ^(x u y jr 0). (3.3) 

The surface density a on Dg- in d3.lt is linearly approximated by 

cr(x,y) ^crij+Sf/x-XiJ+Sljiy-yj), (3.4) 

where (7 !/; = (7(x ;/ y ; ) and 5f- = a x (xi,yj) and 5] j = a y {xi,yj) are constant in the cell Dy. The 

error of the discretization is 0((x— x,) 2 +(y— y ; ) 2 ). Equation d3.4t is the Taylor expansion 
in two dimensions. If a higher order accuracy is required, additional terms in the Taylor 
expansion can be considered. 
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If the surface density is approximated by (13.41) then the force in the x-direction defined by 
03.3|> and (|3.1|) can also be approximated by 
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The evaluation of l|3.5[), (|3.6[) and (|3.7|) can be obtained with the help of the following 
simple integrals, 



( x 2 +y 2)3/2 



dxdy= — ln(y+ J x 2 +y 2 )+C, 



xy 



-dxdy = y\n(x+ \ /x 2 +y 2 ) + C, 



( x 2 + y 2)3/2 

1 



dxdy = — Jx 2 +y 2 + C, 



(x 2 +y 2 ) 3 / 2 



//(x 2 +y 2 )3/2 
The value /C^, • _ ., is equal to 

V% r = ((y-y,) + ^(x-x0 2 +(y-y y ) 2 ) 



rfxrfy = - ^ 2+j/2 +C. 
xy 



' +2 

*•/ 1 
' ~2 



y/4 



(3.11) 
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where the notation g(x) \ =g(b) —g(a). The calculation of /C*l*, :_■< and K^_ v •_ 7 are split 

into two parts by the identity ( x — x\ ) ( x — x v ) = ( x — X{ ) 2 + ( x — x\ ) [x\ — x v ) , and (x — X{ ) (y — 
yy) = (x — Xi){y—yj) + {x — Xi)(yj — yji), respectively. It follows that 



( Xi - x v )lC^ v -_ f + ( (y-ijj) ln(f - Xi + - XiY + (y-ijj) 

toryr)£&j-r+(-y/(*- x iy z +(S-yj) 1 



*•/. i 

1 +2 

X;, 1 



y /4 



* f / l 

' 2 



-V+l 



It is worth noting that the form of F^' , F'j*, and F*' y in (I3.8b - d3.10ll are a type of con- 
volution. It is known that the computational complexity of a convolution of two vectors 
can be reduced to 0(Nlog 2 N) with the help of FFT (see Appendix B). It implies that the 
complexity of this method is 0(N 2 (log 2 N) 2 ). 



4 Self -gravitating force calculation in polar coordinates 

A similar approach is used to develop the method for polar coordinates in this section. 
The singular integral disappears, but the high order of accuracy is not attained because 
there is no explicit form for the integral of an elliptic function. The method in polar 
coordinates is described in detail below. 

The potential function O of fll.lj) under the assumption G = 1 in cylindrical coordinate 
can be expressed as 

&(r,6,z) = -J J J K,{r,r,9,9,z-z)p{r,9,z)fdfd9dz, 
I 

where Kir ,r ,9 ,Q ,z) = — , By dl.41) , the forces on the disk in r- 

K ' ' ' ' ' y / T 2_2rrcos(9-9)+r 2 +z 2 y 

direction and 0-direction become 

— <S>(r,0,O) = / —1C(? r r,9,9,0)a{?,9)rd?d9 (4.1) 

and 

^®(r,9,0) = l J J ^£(f,r,9,e,0)cr(r,9)rdrd9. (4.2) 

We calculate (J4.1|) and (|4.2|) by a numerical approach. 

Since the support of the surface density is compact, contained in a region 1Z = [0,M] x 
[0,2 7t] for some number M > 0, we discretize the radial region in logarithmic form and the 
azimuthal region uniformly as follows. Given a positive integer N, we define A9=2n/ N, 
0</3 <l, /3=/3 o (l-A0), r !+1/2 =/3 N -^ 

& = |(&_ 1/2 +0/+1/2) where ?,/=!,... ,N. It is worth noting that the point r, should be the 
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center of the cell to guarantee the discretization of the surface density is to second order 
and the effect of j8q is to refine the mesh. Here, the proposed method for polar coordinates 
is of first order because a singular integration occurs (see below). Furthermore, the region 
1Z is discretized into the N 2 cells, IZjj = [^-1/2/^+1/2] x [^-1/2/^+1/2] for i,j = l,...,N. For 
j=l,...,N, the cells IZy do not cover the origin and extra cells 7lj= [0,7*1/2] x [^-1/2/^+1/2] 
should be included. For simplification of notation, we denote TZq^ = T^j, 7 = 1, . . . , N. 
The forces in the r-direction and the ^-direction at the point (rj,6j) of the cell TZjj are 

Kr^faW* and F S = ^ M' )- c 4 - 3 ) 

The surface density cr on Rn in (|4.1|) is linearly approximated by 

a{r,d)^a i>j +5l j {r-r i )+5l j {d-e j ), (4.4) 

where = cr(ri,9j) and 5\ ■ = cr r (rj,8j) and d\- = <Jg(rj,8j) are constant in the cell R^j. The 
error of the discretization is 0((r — fj) 2 + (9 — 9j) 2 ). Equation (14.41) is the Taylor expansion 
in two dimensions. 



4.1 The calculation of radial forces 

Let 



r(r, — rcos(0 — OA) 



R <V (f 2 + r 2 - 2fn cos (0 - 9j ) ) ; 



F(r,— fcos(0-0y))(f-r,v) 



R ;y (f 2 + r 2 - 2f r,- cos (0 - ) ) ' 



and 



rfl /" /" f(r i -fcos(9-9 j ))(9-9 jl ) 

/C' e , ..,= // ^ ^dfd0. (4.7) 

' ' ■/ (F 2 +r 2 -2fr,cos(0-0 ; )) 3/2 

The term r, in (|4.6|) is for the formulation of a convolution type. By (|4.1[) and (|4.4|) , we 
have 

N N r r n / \ 

i'=of=i J jR i>,j> or \ j 
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where 



N N 



r(rj — rcos(9 — 9j)) 



F n = EE^';'// '~" v " n " , n drd6 (4.8) 

rfl UK ff f(r,— Fcos(0-0/))(0-0,7) 

Equations (|4.8|) , (|4.9|) , and (|4.10|l can be rewritten as 

N N N 

i'=l;'=l ;'=1 
N N N 

;'=1/'=1 ;'=1 
N N N 

= EE^/;V/+E4'^ r < 4 - 13 > 

l'=l/'=l /'=! 



Let us define F(f,9) = \Jl+f 2 — 2? cos(0), where f is a dimensionless radius. The eval- 
uation of (14.51) , (14.61) and (14.71) can be obtained with the help of the following simple inte- 
grals, 

[ ^~ fCOS W) » = -cos( g )ln(-cos( g ) + ^ + F ( ;0)) + 2CC :^- 1 + C 
J (F 2 +r 2 -2rFcos(0)) 3/2 ^ ' F(|,0) 

:= Hi(-,0) + C 
r 

and 

/ f2(r ~ fCOS(g)) , /2 dr = -rf(3cos 2 (0)-l)ln(-cos(0) + -+F(- / 0)) 



If f 2 f 

+ ^y(-^ cos2 W +3cos W + ^ cos W + ^ ) 

:= rH 2 { r - f 6) + C. 



Following the definition of arid r ;/ we have 

r,- 1+/J n 
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We calculate the value of the integral 

, r0 [ e j'+m fnm f(ri-rcos(6-9j)) 



if ' 



:dfd6 



f -1/2 

7+1/2 



1/2 
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: '+l/2 
: '-l/2 
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■cos(8-0j)ln.(-coB(6-0j)+r/ri+F(r/ri,8-0j)) 

y-1/2 

2cos(0-0 ; )f/r ! -l 

The last integral in the above equation is an elliptic integral and a trapzoidal rule has been 
applied for its evaluation. It is of second order accuracy for the integration of a smooth 
function. Unfortunately, the presence of a singular function in terms of ln(l — cos(0)) 
reduces the accuracy of the proposed method for polar coordinate to first order. Finally, 
the value /Cjf jV • is approximated as follows and is used in the numerical calculation, 

1 

~ 2 C^l( r i'+l/2/ r i'fy+l/2— Qj) -Hi(r j >_ 1 / 2 /ri,6ji +1 /2-9j) 
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where the notation /(•)]„ = |(/(fc)+/(fl))(&— fl). Similarly, 
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4.2 The calculation of azimuthal forces 

-0,0 ir6,r 



Next, we introduce the calculation for JC u i AJ f K?S_ ., • _•„ and /C^f f/ ■_ ., 
In particular, we calculate the value of the integral 
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and 



v6,6 

K i-i',i-f 



-i+cos 2 (0-e y ; 



sin(0 - 6j) ( 2cos 2 (6 - y ) - +cos(0 -6,)) 
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5 Order of accuracy and a comparison study 
5.1 Order of accuracy 

We investigate the numerical errors induced by the truncation introduced in ( 13.41 ), which 



is 



X — Xj 

V/ ((i-x^+^-y,)*) 4 



o(((Ax) 2 + (Ay) 2 )//^ -I ■ " - rfafr . 



The last integral in the above estimation is 0(|lnAx|) which gives the numerical error 
of order 0((Ax) 2 |lnAx|) = 0((Ax) 2 ~) with Ax = Ay. Three types of norm are used to 
measure the errors between the numerical and analytic solutions. The L 1 , L 2 , and L°° 
norms of a function / on a domain O are defined as 

r V /p 

J^\f(x)\ p dxj , for p = l,2, and ||/|| 00 = esssup n |/(a:)|. 

The errors between the analytic and numerical solutions for various resolutions using 
different norms L 1 , L 2 , and L°° demonstrate the convergence in total variation, energy 
and pointwise senses, respectively. 

We verify that the proposed method is of second order accuracy by demonstrating the 
following examples, a D2 disk Q, a non-axisymmetric disk consisting of two superposed 
D2 disks and a non-axisymmetric disk describing a pair of spirals. 

Example 1 . The D2 disk has the surface density 

EnfM-l ^(l-tf/" 2 ) 3/2 iorR<a, 



where R = y x 2 +y 2 and a is a given constant. The corresponding potential on the z = 
plane is 



O D2 (R,0;a) 



^^( 8a 4_ 8a 2 R 2 +3R 4) iorR<a 

for R > cl, 
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Table 1: This table demonstrates the errors and order accuracy of the proposed method for the D2 disk for 
various number of zones N = 2 k from k = 5 to 10. It shows that the order for the D2 disk is about 1.8 or 1.9 
order in L 1 and L 2 norm. 



and the radial force is found as 



I 8a 3 



^^(Acc 2 -3R 2 ) for R<« 



R(4a 2 -3R 2 )sin- 1 (f)-a(2a 2 -3R 2 ) v / l-a 2 /K 2 



for R > a. 



Without loss of generality for studying the order of accuracy let us set the computational 
domain Q= [— 1,1] x [— 1,1], <7n = G = 1 and a = 0.25. We illustrate the contour plots of the 
surface density, x-directional force, y-directional force, radial force, residuals between 
analytic and numerical solutions for x , and radial directions for N = 1024 in Fig. [TJ The 
residuals show that the largest errors occur in regions surrounding the edge of the disk 
where the second derivative of the surface density with respect to radius is infinite. 

In Table [TJ the column E% and is the error of the x directional force and R radial 
direction by p-norm, p = 1,2, and 00, between the analytic and numerical solutions. The 
column Oj in Table [JJis the order of accuracy as measured by log 2 (Ex (2 k ~ 1 )/E x 1 (2 k )) for 
k = 6 to 10 and similarly for Oj^. These errors show that this method is nearly of second 
order accuracy. More precisely, we obtain the order of convergence to be about 1.8 or 1.9 
as measured by the L 1 and L 2 norms for the simulation of a D2 disk. The L°° norm only 
demonstrates the convergence, since the second derivatives of the surface density of the 
D2 disk are not bounded. 



We continue to use the D2 disk as an example and a unit disk D (0,1) =Q= [0,1] x [0,27r] 
as the computational domain to investigate the self-gravitational force in polar coordi- 
nates. The value /3n = 0.99 is set. We show the contour plots of the surface density, radial 
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-0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 
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Figure 1: The numerical solutions of a D2 disk for N = 1024, the contour plots are surface density (upper 
left), the y-directional force (upper right), the x-directional force (middle left), the radial force (middle right), 
the difference between analytic and numerical solutions in x direction (lower left), and the difference in radial 
direction (lower right). The values in the lower contour plots are the absolute difference in the common 
logarithmic scale. 
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1.002 
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3.646E-2 
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1.035 


1.000 


256 


1.754E-2 
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3.403E-2 


128/256 


1.056 


0.947 


1.000 


512 


8.762E-3 


1.049E-2 


1.701E-2 


256/512 


1.001 


1.000 


1.000 



Table 2: This table demonstrates the errors and order accuracy of the proposed method for the D2 disk for 
various number of zones N = 2 k from k = 5 to 10 on polar coordinates. It shows that the order for the D2 disk 
is about 1 in each norm. 



Surlace Density a(x.y) 




-0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 -0.8 -0.0 -0.4 -0.2 0.2 0.4 0.6 0.8 -0.8 -0.0 -0.4 -0.2 0.2 0.4 0.6 0.8 



Figure 2: The numerical solutions of a D2 disk for N = 512 to investigate the self-gravitational force calculation 
in polar coordinate. From left to right, the contour plots are surface density, the radial directional force, the 
difference between analytic and numerical solutions, respectively. The values in the right contour plot are the 
absolute difference in the common logarithmic scale. 

force, and the difference between analytic and numerical solutions for N = 512 in Fig. [2] 
and the order of accuracy is only about 1 as given in Table |H The largest errors occur in 
regions not only surrounding the edge of the disk, but also close to the origin. Although 
the surface density at the origin is smooth, the singular elliptic integral introduces signif- 
icant error there. Hereafter, we concentrate on the self-gravitational forces in Cartesian 
coordinates. 

Example 2. The disk 02,2 of two superposed D2 has the surface density £d 22 =2 ^d 2 + 
^d 2 (^2;«)/ where R\ = y 7 (x — l/4) 2 +y 2 and R2 = \J (x + l/4) 2 +y 2 . This example repre- 
sents a non-symmetric distribution of the surface density of a disk. The results are shown 
in Table |3] and Figure |3] This result is similarly to Example 1. The factors 0~ of errors 
in Table [3] are non-mono tonic as the numerical resolution, N^, increases. This may be 
due to the distribution of the surface density on grid cells, the centers of which can shift 
with varying numerical resolution. However, the total variation and energy shows the 
convergence and the order of accuracy is about 1.8 and 1.9 respectively. 
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-0.8 -0.B -0.4 -0.2 0.2 0.4 0.6 0. 



Figure 3: The numerical solutions of a D22 disk for N = 1024, The top contour plot is the surface density. 
The contour plots in the second row are the x-directional, y-directional, and radial forces, respectively. The 
corresponding errors between the numerical and analytic solutions in the third row. The values in the contour 
plots in the third row are the absolute errors in the common logarithmic scale. 
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Table 3: This table demonstrates the errors and order of accuracy of the proposed method for the D2 2 disk for 
various number of zones N = 2 k from k = 5 to 10. It shows that the order for the disk is about 1.8 or 1.9 
order in L 1 and L 2 norm. 

Example 3. As another example of a non-axisymmetric potential, we consider a logarith- 
mic spiral disk. Since an analytic pair for the surface density and potential are not known, 
we assume a surface density profile of the form 

L LS (r,0)=e- 2r2 (2+cos(20 + 16r)). 

To investigate the order of accuracy, the solution at the finest mesh size is regarded as 
the true solution. For various coarser resolutions, the values at some specific position are 
taken to be the average of the four closest to the position. The results are shown for the 
method based on Cartesian coordinates in Table H] and Figure HJ It can be seen that the 
order of accuracy is about 1.5 for the L 1 norm and about 1 for the L 2 norm. The L°° norm 
is only convergent. 

5.2 A comparison study 

The goal of this paper is to calculate the self-gravitational forces with as few restrictions 
as possible. The most straight forward approach is to solve for the potential via dl.lD 
and obtain the self-gravitational forces by taking its derivatives. If one uses the finite 
difference or finite element method on dl -lb , the discretization is 

- < t > i+l,j,k+2®i,j,k- < S > i-l,j,k -^i,i+l,k + ^i,j,k-^i,j-U -®i,jMl+ 2 ®i,j*-®i,j*-l _ _r 
(Ax) 2 + (Ay) 2 + (Az) 2 ~ k 

where O, ^ / c =0(x;,y ; -,zjt) and fi,j,k=f{ x i'yj, z k) based on the uniform mesh grids (j;,y ; ,z/ c ). 
Here, /,• = for k 7^ 0. For such an approach, artificial boundary conditions should be 




Figure 4: The numerical solutions of a logarithmic spiral disk for N = 512 to investigate the self-gravitational 
force calculation. The contour plots illustrate the surface density (upper left), x-force (upper right), y-force 
(low left), and radial force (lower right). 
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Table 4: This table demonstrates the errors and order of accuracy of the proposed method for the spiral disk 
for various number of zones N = 2 k from k = 5 to 9. It shows that the order for the spiral disk is about 1.5 or 
1.0 in the L 1 and L 2 norms, respectively. 

imposed and a fully 3-dimensional calculation must be undertaken. We point out that the 
l ll.lt can not be reduced to the two dimensional numerical partial differential problem, 
i.e., 

-®i+i,jjo+2®i,jfl-®i-i,jfl -®i,j+ifl+2®i,jfl-®i,j-ifl _ , 

(Ax) 2 + (Ay) 2 " W 

Any numerical solution of the partial differential problem will involve 0(N 3 ) unknowns. 
It follows that the linear complexity of such an approach, viz. multigrid method, is at 
least 0(N 3 ). For an infinitesimally thin gaseous disk problem, this approach does not 
appear to be suitable. 

Alternatively, one can solve the reduced equation given by 

*(z,y,0) = -G / / dxdy 
J J ^/(x-x) 2 + (y-y) 2 

or 

®(r,e,0) = -G ( ( a{j '^ TdrdO. 

J J ^r 2 +r 2 -2rrcos(6-6) 

In this case, one can consider using bases functions on a two dimensional space as in 
a spectral method. Unfortunately, this approach requires a treatment for the boundary 
conditions. A possible way to deal with this issue is to impose periodic boundary con- 
ditions. However, it is not realistic for a gravitational force calculation because gravity 
is a long range force and not periodic. As an alternative, a method without the periodic 
assumption has been proposed for polar coordinates PJ. The approach in (3] transforms 
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the polar coordinate (r,0) into the coordinate (u,8) by setting r = e" or u = \n(r). The 
potential-density pair in term of the reduced surface density and reduced potential is 
given in (3], and it is 

1 f°° 



and 



1 r°° 

e ll/2 <P(e u ,6) = - — G'£ j K(ot,m)A m (cc)exp[i(me+ctu)]dcc, (5.2) 

ATI J —oo 



where K is real and positive and is defined as 

, _ 1 r[(m + l/2+k)/2]r[(m + l/2-m)/2] 
[a,m) ~ 2 r[(m+3/2+k)/2]r[(m+3/2-m)/2] ' 

We regard this method as one of the spectral methods because Fourier series e~"" 6 and 
Fourier integral e~ mu are used. To apply this method to the D2 disk using the polar 
coordinates, we transform the bounded unit disk D(0,1) = [0,1] x [0,2 7r] to the unbounded 
domain U — (— oo,0] x [0,2 7r]. In this special case, we only need to compute m = and 
truncate 

A (ct)= f° e 3u/2 cr(e u )e- iau dun f° e 3u/2 a(e")e- iau du, (5.3) 

J -co Ju min 

where the value u m i n is to approximate —00. The truncation produces a hole in the unit 
disk and can introduce significant errors at the origin. Given a positive integer N and 
base on the discretization for the radial region in the previous subsection, to calculate 
(J5.3|) and (|5.2|) by the trapzoidal rule. The variation of the potential with respect to ra- 
dius is illustrated in Figure |5] The profile on the left panel shows that the numerical 
and analytic solutions for the Kalnajs' method agree well except close to the origin for 
N — 1024. The small window embedded within the panel zooms in on the residuals be- 
tween numerical and analytic solution on the interval [0,0.3]. It is seen that the truncated 
portion contributes to significant errors near the origin. In contrast, the application of our 
proposed method to the calculation of potentials leads to the results shown in the right 
panel of Figure Although the singular integration still remains due to the unbounded 
domain, our proposed method on either Cartesian and polar coordinates is preferable 
since a hole near the origin is not introduced. 

Finally, a third approach is to directly calculate the integrals and obtain the potential. 
For any given mesh grid, the total amount of complexity is 0(N 4 ) based on the num- 
ber 0(N 2 ) of mesh zones. If we restrict ourselves to a uniform grid and use the FFT 
technique, the complexity can be reduced from 0(N 4 ) to 0(N 2 ). In other words, a fast 
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Figure 5: The variation of the potential with respect to radius using Kalnajs' method (left) and the proposed 
method (right). The residuals are shown in the small window in each panel and show that the Kalnajs' method 
have significant errors near the origin, which are eliminated in the proposed method. 

algorithm of linear complexity is obtained. It is common to start with 

O(x,i/,0) = — G J J lC(x — x,y—y,0)a(x,y)dxdy 

N N , , 

= -G£ / / K(x-x,y-y,Q)a(x,y)dxdy. 

i=li=l J JD hi 



and to introduce a softening parameter e to approximate 

G 



I Id 1C ^~ x, y~ y ^ 



a{x,y) ; 



^e 2 + (x i ,-x i ) 2 + (y r -y j y- 



■. / / a(x,y)dxdy. 



Since the goal is to calculate the forces, the order of accuracy is reduced when taking 
the numerical differentiation on the numerical solution of potentials. For polar coordi- 
nates HI, the value of JC is approximated by 



2 (cosh (iiji — u {) — cos (Qji — 6j ) ) 



where M,/ =ln(x ! v) and Uj = ln(x ; ). Note that when (i',f) = (/,/), /C is undefined. An 
approach to avoid the singularity problem can be found in (TJ. On the other hand, the 
proposed method avoids the singularity problem by directly evaluating the forces, hence, 
raising the order of accuracy. For Cartesian coordinates, we choose the softening param- 
eters as the mesh size e = Ax. The errors for the disks D2 and Di,i are shown in Table |5] 
and Tabled respectively. It reveals that the accuracy when using the softening parameter 
approach for the D2 and D22 disks is of first order in the L 1 and L 2 norms. For the L°° 
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Table 5: This table demonstrates the errors and order accuracy of the softening parameter method for the D2 
disk for various number of zones N = 2 k from k = 5 to 10 in Cartesian coordinates. It shows that the accuracy 
for the T>2 disk is about first order. 



norm, the order of accuracy for the D2 disk is about 1. For the 02,2 disk, this method loses 
accuracy. In comparison with our proposed method for Example 1 and Example 2, our 
methods are more accurate and the order of accuracy is verified. 

We implement the proposed method using MATLAB 7 software under the computer 
system, Intel Core 2 Duo CPU 1.8GHz with 2 GB RAM. The CPU time measurement 
information of the proposed method is compared with the direct method in Table [7] We 
list the CPU times in evaluating the kernels /C ' , the force calculations of convolutions, 
and the whole process. The measurement is evaluated by the mean of 40 simulations. It 
shows that the CPU times of both of the proposed method (P.M.) and the direct method 
(D.M.) are comparable. 



6 Discussion and conclusion 



We have presented a near second order method for calculating the self-gravitating force 
of an infinitesimally thin disk for Cartesian coordinates. For polar coordinates, we find 
that the method is near first order, ~ 0.89, only. To quantify the accuracy, we define 



rh 1 

/ ln(l-cos(0))d0--(ln(l-cos(0 J t))+ln(l-cos(-0 J t))2(0 Jt , 

J-e k 



where 9^ = 1 /2 k . Table [8] reveals that the accuracy of the trapzoidal rule for the integration 
of the function ln(l — cos(0)) is nearly of first order. With an improvement of the singular 
integration of ln(l — cos(#)), the accuracy can be increased for the proposed method in 
polar coordinates. 

We note that the fast Fourier transform is only used to reduce the computational time. 
For the practical computation, one can extend the range of the summation in (|3.8|) . By 
setting = whenever either i' or f is in the range — N+l to 0, the value of any of 
the F*j° is unaffected. Furthermore, we can take cfy^/ to be periodic since the sum 03.8[) 
does not involve any values of i 1 and / outside the first period. We are also free to take 
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Table 6: This table demonstrates the errors and order accuracy of the softening parameter approach for the 

D2,2 disk for various number of zones N = 2 k from k = 5 to 10. It shows that the order for the Diz disk is about 

first order in L 1 and L 2 norm. For measurement of L 00 norm, this method may fail in convergence under the 
pointwise sense. 
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Table 7: This table demonstrates the CPU time measurement of the proposed method (P.M.) and direct method 
(D.M.) with softening parameters. The whole process consists of the generation of kernels and the forces of 
calculations. It shows that the CPU times of both of P.M. and D.M. are comparable. 
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Table 8: This table demonstrates the accuracy of the trapzoidal rule for the integration of the function ln(l — 
cos(8)) is near of first order ~0.89. 
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/C^ljv • v periodic by defining it to be the periodic function that agrees with \6.7) for i — i' 
and j—f in the range [— N+1,N] of the Green function. 

An important feature of our approach is that the boundary is not assumed to be pe- 
riodic. Our approach is limited to the Cartesian and polar coordinates with uniform and 
logarithmic grid discretization, respectively which allows for rapid computation. That 
is, the restriction of a convolution of two vectors provides the rapid computation, but it 
is restricted to a grid discretization that is either uniform or logarithmic. If the discretiza- 
tion is arbitrary, then the FFT is not suitable. 

We point out that our method may be useful for gravity computations on a nested grid 
consisting of uniform grids having different grid spacing designed to resolve a central 
region with a finer grid. Such an approach would be complementary to the fast algorithm 
for solving the Poisson equation on a nested grid presented by Matsumoto and Hanawa 
0. 
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Appendix A: The calculation of the force in the y-direction in 
Cartesian coordinate 



Let 



and 



By d3D and ffll) . we have 
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m,i-r= L 91 L.n dxdy- (6-3) 
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The evaluation of l|6.1[), (|6.2[) and (|6.3|) can be obtained with the help of the following 
simple integrals, 
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The value /C^ -, ■-, is equal to 



K\_ ViH , = - In ( (x - Xi ) + ^x-Xif + O-yjY) 



*•/ i 
' ~2 



(6.7) 



The calculation of Kf^ v and K?^, are split into two parts by the identity {y—yj) (y- 

y/')=(y-y;) 2 +(y-y/)(y/-y;')/and (y-yjO(*-^)Hy-y;0(*-^0+(y-y/)(*i-^')/« 

spectively. It follows that 

= {yj-y?)1C\% iHI + ( (x - X,-) ln(y- y y + ^/(y - y ; -) 2 + (f - X/) 2 )) 
= - x v )Kf v ._ f + (- ^ (y-yy) 2 + (x - x,) 2 ) 
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Appendix B: Calculations of convolution of two vectors by FFT 

It is known that the FFT of a vector u n , n = —N,...,N—1 can be rewritten as 

N-l 

u k = £ u n e-> 27Tkn/2N ,fork=-N,...,N-l, 

n=-N 

and its inverse FFT is given by 

-i N-l 

M « = ^j E V 27r " n/2N ,forn = -N,... / N-l. 

/iV jc=-N 

Let us consider two vectors u n , n = 0,...,N—l and v„, n = —N+l,...,N—l and their inner 
product 

N-l 

w k = u nVk-n, for fc=0,...,N-l. 

»=0 



We set zc?fc = 0, fc= —N,...,0 and 



N-l N-l N-l 

fc=-N k=-Nn=-N 

N-l N-l 

= ^ u n e~i 2nmn/2N J^j v k n e~i 2nm ( k ~ n ^ 1N 

n=-N k=-N 
N-l N-l 

= £ Mn£ >-^ m,3/2N X] i> fc <T' 27rmfc/2N 

k=-N fc=-A7 
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This gives us 

w m = u m -v m , for m = —N,...,N— 1. 

Applying the inverse FFT on the above equation, we recover the vector zv m , m=—N,...,N— 
1. The vector «;,„, m = 0,...,N— 1 is the desired result. 



